/*
 * InlineStream.cpp
 *
 *  Created on: Jan 5, 2011
 *      Author: blyth
 */

#include "InlineStream.h"
#include "PeakReference.h"
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/RandGauss.h"
#include <cmath>
#include <iostream>
using namespace std;

InlineStream::InlineStream(unsigned int length, unsigned int numPeaks, long seed, double* refPeak,
		unsigned short refSize) :
	length(length), numPeaks(numPeaks), stream(new unsigned short[length]), refSize(refSize), refPeak(refPeak) {

	randEngine = new CLHEP::MTwistEngine(seed);
	randGauss = new CLHEP::RandGauss(randEngine);
	peakInfo = 0;

	reRandomize();
}

InlineStream::~InlineStream() {
	if (stream)
		delete stream;
	if (peakInfo) {
		for (int i = 0; i < numPeaks; i++)
			delete peakInfo[i];
		delete peakInfo;
	}
	delete randGauss;
}

void InlineStream::reRandomize() {
	if (peakInfo) {
		for (int i = 0; i < numPeaks; i++)
			delete peakInfo[i];
		delete peakInfo;
	}

	for (int i = 0; i < length; i++)
		stream[i] = 500;

	peakInfo = new double*[numPeaks];
	for (int i = 0; i < numPeaks; i++) {
		peakInfo[i] = new double[2];
		peakInfo[i][0] = randEngine->flat() * 15884;
		peakInfo[i][1] = randEngine->flat() * (length - refSize);
	}

	addNoise();

	PeakReference PR(refPeak, refSize);
	for (int i = 0; i < numPeaks; i++) {
		int start;
		double offset = peakInfo[i][1] - (int) peakInfo[i][1];
		if (offset == 0)
			start = (int) peakInfo[i][1];
		else
			start = (int) peakInfo[i][1] + 1;
		for (int j = 0; j < refSize && start + j < length; j++)
			stream[start + j] += PR.cubicInterpolationAt(j + offset) * peakInfo[i][0];
	}

}

void InlineStream::addNoise() {
	for (int i = 0; i < length; i++)
		stream[i] += 500 + (int) (randGauss->fire(0., 2.65));
}
